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Abstract 



Fixed points are fundamental states in any dynamical system. In the case of gene regulatory networks (GRNs) 
they correspond to stable genes profiles associated to the various cell types. We use Kauffman's approach to model 
GRNs with random Boolean networks (RBNs). We start this paper by proving that, if we fix the values of the 
source nodes (nodes with in-degree 0), the expected number of fixed points of any RBN is one (independently of 
the topology we choose). For finding such fixed points we use the a-asynchronous dynamics (where every node 
is updated independently with probability < a < 1). In fact, it is well-known that asynchrony avoids the cycle 
attractors into which parallel dynamics tends to fall. We perform simulations and we show the remarkable property 
that, if for a given RBN with scale-free topology and a-asynchronous dynamics an initial configuration reaches a 
fixed point, then every configuration also reaches a fixed point. By contrast, in the parallel regime, the percentage of 
initial configurations reaching a fixed point (for the same networks) is dramatically smaller. We contrast the results 
of the simulations on scale-free networks with the classical Erdos-Renyi model of random networks. Everything 
I indicates that scale-free networks are extremely robust. Finally, we study the mean and maximum time/work needed 

^. to reach a fixed point when starting from randomly chosen initial configurations. 

in 
On 
in 1 Introduction 

A Random Boolean Network (RBN) is a model of a GRN. This model was introduced by Kauffman in 1969 ll26l and 
it corresponds to a directed graph composed by N genes (nodes) where each of these genes can be either expressed 
(state 1) or not expressed (state 0). Each gene v receives K randomly chosen genes as input (K nodes pointing to 
v). In other words, the dynamics is defined locally. The idea of Kauffman was to choose, independently for each 
node v, any of the 2 2 possible functions f v : {0, 1} K — > {0, 1} with equal probability. Finally, a global dynamics 
/ : {0, 1}^ — > {0, 1}^ was defined by applying the local functions in parallel (sometimes this is referred to as 
synchronous updates). Kauffman found, through simulations, the existence of a phase transition at K = K c = 2 
(from an ordered phase to a chaotic phase). He also suggested that the number of attractors grows exponentially with 
N in the chaotic phase, proportional to N in the ordered phase, and proportional to y/N at the critical value K c = 2. 

Kauffman's model (and some natural generalizations) motivated a lot of work carried out mainly by physicists. 
They tackled problems arising from Kauffman's definitions both analytically and numerically: the critical in-degree 
value f8l [T6ll . the number of attractors lfTTl[T6 l, the length distribution of the attractors [4 10], etc. 

From the biologists point of view a fundamental challenge of the post-genomic era is the possibility of simulating 
the dynamics of real genetic networks. The enormous amount of available data of molecular interactions within 
the cell made it possible to examine critically the original RBN model. The first observation was that the topology 
assumption was inadequate. Contrasting the uniform topology assumed in the original RBNs, it has been shown that 
real genetic networks exhibit a scale-free topology |3] [12] [20] [28) . In such topologies a small fraction of the genes 
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are highly connected whereas the majority of the genes are poorly connected (HBO). Therefore, RBNs dynamics with 
scale-free topology started to be intensively studied ll5l [131 24 1 . 

Another criticism towards Kauffman's model was that nodes were updated in parallel. Experimental results con- 
firmed a rather intuitive fact: that genes transition between expressed and non-expressed states at different times ifTTl 
WD . Informally, there is no global clock that allows transition to happen only at ticks. Therefore, RBNs with scale- 
free topology and asynchronous dynamics is a natural model to be analyzed 1 14 1. We would also like to point out that 
asynchrony in the classical RBN model has also been studied in [21 22, 27] and l29l . 

Our main contributions. Fixed points are fundamental states in any dynamical system. In the case of GRNs 
they correspond to stable gene expression profiles associated to the various cell types. This interpretation has been 
used for modeling europhil differentiation |23|, expression patterns of the segment polarity genes in Drosophila 
melanogaster [2|, flower organ specification in Arabidopsis thaliana (6], etc. 

The goal of this paper is to study the of existence fixed points in RBNs with scale-free topology. Notice that a fixed 
point does not depend on the timing (synchronous/asynchronous) of the update rule. In fact, in a fixed point every 
gene is in a stable state with respect to its local input. Therefore, there is no way to change such global configuration. 
In that sense, a fixed point is a very robust object of a network, despite the fact that its basin of attraction can change 
depending on how the update rule is implemented. 

Once the topology of the network is (randomly) generated some nodes will have only outgoing arcs. These were 
called source nodes by Albert [3| and we will use the same terminology here. Since their states do not depend on the 
state of any other node we fix their values arbitrarily. In this paper we prove that, for every given assignment to the 
input nodes, the expected number of fixed points is exactly one (for every topology). 

Now the main (and natural) question arises. Given a RBN, how can we actually find its fixed points? It is clear that 
testing all the 2 N configurations is impossible. For answering the question we come back to the issue of asynchrony. 
In fact, in 1994 Bersini and Detours studied an asynchronous version of the cellular automaton Game-of-Life |9|. 
They observed that the introduction of asynchrony modified the dynamics from a behavior with long transients to a 
behavior with fixed points. This is rather intuitive: asynchrony is a way to avoid the cycle attractors the deterministic 
(parallel) implementation tend to fall into. 

Therefore, the strategy is clear. Given a RBN, we implement the a-asynchronous dynamics for different values of 
ot (0 < a < 1) [19]. Roughly, this means that, at each time step, each gene is updated independently with probability 
a. When a varies from 1 down to the dynamics evolves from the fully deterministic synchronous regime to a more 
asynchronous regime. When a = 0we choose randomly only one node at each step. 

Our simulations show that RBNs with scale-free topology for which there exist a fixed points every initial con- 
figuration converges to a fixed point when < a < 1. On the other hand, when a = 1, the percentage of initial 
configurations that reached a fixed point varies greatly form one network to another. In some cases the percentage is 
close to 0. In average (considering all the networks we use) the percentage is ^28.9%. 

In order to find properties which could be associated exclusively to the topology of the network we compare 
the results of the simulations on scale-free networks with the results on the classical Erdos-Renyi model of random 
networks [18]. The main difference with the scale-free topology is that here we generated networks for which the 
percentage of initial configurations converging to a fixed point is close to 0. Finally, we study the mean and maximum 
time/work needed to reach a fixed point when starting from randomly chosen initial configurations. 

Therefore, a remarkable and distinguishable dynamical property arise on RBNs with scale-free topology: robust- 
ness of convergence under asynchronous update. This fact could provide some insight about why such topologies are 
ubiquitous in GRN. Furthermore, asynchronous updating could be a natural mechanism present in GRNs in order to 
avoid cyclic dynamics [ 29 1 . 



2 Network model 

A Random Boolean Network (RBN) corresponds to a directed graph composed by N genes (nodes) where each of these 
genes can be either expressed (state 1) or not expressed (state 0). We will refer to the nodes of a RBN as Ui,V2, ■ • • ,Vn. 
We define K\ n as the in-degree of node Vi and K° ut as the out-degree of node Uj. Every zero in-degree node is called 



a source node, while every non-zero in-degree node is called an internal node. A configuration of the network is a 
vector s € {0, 1} N that associates a binary state to each of the nodes. 

2.1 Dynamics 

We assign to each gene Vi a local transition rule fa : {0, l} K i +1 — *■ {0, 1}. Informally, the value of fa depends 
on the state of the K\ n input nodes together with the state of Vi itself, source nodes always remain in the same 
state. More precisely, if K\ n = 0, then fa(0) = and fa(l) = 1. For each internal node Vi we construct ran- 
domly its local transition function as follows. Call k the in-degree of Vi (i.e, k — K\ n ). There are 2 2 possible 
functions of the form / : {0, l} fe — ► {0, 1}. A straightforward approach is to choose one of these functions from 
a uniform probability distribution. Nevertheless, before selecting a function, we should rule out those which do not 
strictly depend on all of its arguments (otherwise we would not be respecting the network topology). To define this 
concept precisely, we say that a function / : {0, l} fe —J- {0, 1} strictly depends on its arguments iff for all j <G 
{1, 2, • • • , k}, there exist a;i,X2, ■ ■ ■ , Xj-i, Xj+i, ■ ■ ■ , Xk € {0, 1} such that f(x\, Xi, • • • , Xj-i, 0, Xj+i, • • • , Xk) =/= 
f(xi,X2, ■ ■ ■ ,Xj-\, l,Xj+i, ■ ■ ■ ,Xk)- We fix function fa by selecting one (randomly) among all those functions 
strictly depending on all of its arguments. 

The global update rule is characterized by a real parameter a <E [0, 1]. We denote the global transition rule by 
<l? a : {0, 1} N — *■ {0, 1}^, and we define it through the following protocol: 

1 . Select each internal node independently with probability a. We call selected nodes this set of randomly chosen 
internal nodes. For the special case a = we select randomly a single internal node (i.e, the set of selected 
nodes is a singleton). 

2. Update in parallel all the selected nodes (applying the local transition rule in all the nodes belonging to such 
set). Do not change the state of the other nodes (input nodes and non-selected nodes). 

Let s t £ {0, 1}^ be a configuration of a RBN at time t € N. A stochastic trajectory, starting from the initial 
configuration Sq, is the sequence so, Si, S2, . . ., where Sj = $ a (sj_i). 

The parameter a can be thought of as a measure of parallelism in the update process. The strictly sequential- 
random policy rule is captured by a = 0, where only one internal node is updated at each step. Similarly, by making 
a = 1, we represent the full parallel-deterministic policy rule, where all the internal nodes are updated at each step. If 
a = 1 we call the resulting trajectory the deterministic trajectory of the system. A configuration s is a fixed point of a 
RBN iff Pr{<J> Q (s) = s}=l. This is equivalent to say that $i(s) = s. Therefore, s is a fixed point regardless of the 
choice of a. 

We measure time simply by counting the applications of <J> Q . Therefore, the time to go from So to sy in a trajectory 
is T. This notion of time neglects the fact that the computational effort to evaluate $ Q depends on a. The expected 
number of fas to be evaluated is aNj where Nj is the number of internal nodes. Thus, we define the work to go from 
state s to s<r to be aT for all a > 0. For the special case a = 0we define the work as jj-T. 

When considering the dynamics of the network with a = 1, it is clear that the trajectory will eventually become 
cyclic. Note that if the system reaches a fixed point, the length of the period is 1. Both the time to enter the cycle and 
the period are bounded by 2 Nl , where iVj is the number of internal nodes. However, if a < 1, the idea of cycle is 
insufficient to describe a trajectory that does not reach a fixed point. Instead, we have a set of configurations which are 
visited infinitely often. This greatly complicates the numerical determination of any eventual convergence to a fixed 
point. We therefore select a somewhat arbitrary time horizon to interrupt our simulations. 

Notice that, given an initial configuration so, the set of configurations reachable from So by repeated applications 
of $i is a subset of the set of configurations potentially reachable from So by repeated applications of <fr Q when a < 1. 
In other words, the stochastic trajectory can visit a larger portion of the state space than the deterministic one. 

A remarkable feature of selecting the local update functions randomly, as we did, is that we can predict, in a 
statistical sense, the number of fixed points in the network. Consider an arbitrary configuration s e {0, 1}^. What is 
the probability of s to be a fixed-point? First notice that s is a fixed point if and only ifsi = fa(si 1 ,...,Si k .) for every 
i E {!,..., N}, where v^, . . . , Vi k . are the input nodes of i>,. 



It follows from the definition of the (pi's functions, that Pr{0j(sj i; . . . ,Si k ) = 0} = Pr{^i(si 1 , . . . ,Sj fc ) = 1} = h 
for every internal node Uj. The reason for that is the following: if a function (pi strictly depends on all of its arguments 
then the complementary function (the one where we replace l's with 0's and 0's with Is) also does. On the other hand, 
Pr{(pi(si) = s^ = 1 for every source node Vi. Therefore, the probability of s to be a fixed-point is 2~ Nl , where Nj 
is the number of internal nodes. 

Let X s be the Boolean random variable that equals 1 if the configuration s is a fixed point and otherwise. The 
expected number of fixed points is 
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Note that N — Nj = Ne is the total number of source nodes of the network. If we fix all these source nodes to 
some arbitrary value in the initial configuration, then the expected number of fixed points goes down to 1 . In fact, in 
that case, the number of effective different configurations is 2 n ~ Ne = 2 Nl and the expected number of fixed points 
becomes 2 Nl ~ Nl = 1 instead of the original 2 N ~ Nl . 

Until this point, neither the definitions nor the theorem we just proved make assumptions about the topology of the 
RBN. In the following subsection we will describe two families of topologies that have been studied in the literature. 

2.2 Topology 

We start describing here a process by which we construct a directed scale-free network with N nodes and average 
in/out-degree equal to k. Let iVo and k be positive integers such that k < Nq < N. The process starts from a directed 
clique with No nodes (i.e, Nq(Nq — 1) arcs). Call these nodes vi,vz,..., vn . The process now involves N — No 
growth stages, numbered No + 1, No + 2, . . . , N. At each stage, a single node is added to the network. 

Call Vi the node added at stage i. We will also add k edges to the growing network. We toss a fair coin and proceed 
as follows: 

1. In the case of heads we add k edges pointing from k different nodes in {vi, . . . Vj-i} towards Vi. These k nodes 
are selected randomly following a preferential attachment rule such that the probability of Vj to be selected is 
proportional to K° ut + 1. 

2. In the case of tails we add k edges pointing from Vi to k different nodes in {u 1; . . . Uj_i}. These k nodes are 
selected randomly following a preferential attachment rule, such that the probability of Vj to be selected is 
proportional to K" 1 + 1. 

We will refer to the process just described as the BA algorithm because it is based on previous work by Barabasi 
and Albert [7|. Following 11241 . we started with a clique of size No — 5 and average in/out-degree k — 2 to create 
the topology using the BA method. If a RBN has a topology created by the BA algorithm, we will call it a scale-free 
network. 

The BA method, as defined here, never creates an edge from some node to itself. The same assumption is present 
in Q and [24] and is pervasive in the literature. It is noteworthy that the theorem proved in Subsection 2.1 does not 



require a special topology. Therefore, it still applies to networks where auto-regulation is present. Thus, for given 
values for the input nodes, the expected number of fixed points would still be 1 . It is conceivable that the probability 
distribution of the number of fixed points will change, though. 

In the next sections, we compare the results on scale-free networks with the classical Erdos-Renyi model of ran- 
dom networks [18] (random directed graphs). By doing so we intend to find properties which could be associated 
exclusively to the topology of the network. Random networks of this well-known Erdos-Renyi family can be easily 
generated. Let p be a real number in the [0, 1] interval. Each potential link (vi, vf) is selected independently with 
probability p. Therefore, the expected number of links in the Erdos-Renyi network is pn(n — 1), where n is the num- 
ber of nodes in the Erdos-Renyi network. We will call this process the ER algorithm. If a RBN has a topology created 
by the ER algorithm, we will call it an Erdos-Renyi network. 

Conceptually, to make such comparison "fair," we have to adjust the parameters used by the generation algorithm 
so the networks generated have some common statistical properties. It seems reasonable to preserve both the number 



of internal nodes and the average in-degree. The first parameter determines the number of states the system can 
be in, after the state of the source nodes have been fixed. The average in-degree is a measure of connectivity and 
density. Formally, if we run the topology generation BA algorithm presented for scale free networks with parameter 
N (recall that N = 5 and k = 2) we will obtain a scale-free network Gba- The average in-degree is k — 2. Call 
Ni the expected number of internal nodes of Gba- Similarly, if we run the topology generation ER algorithm with 
parameters n and p, we will obtain a graph Ger- 

The problem we wish to solve is: Given Ni find n and p such that the expected in-degree of nodes in Ger is 2 
and the expected number of internal nodes of Ger is Nj. Therefore, 

p-(n-l) = 2 (1) 

n ■ (1 - (1 - p) n - 1 ) = Nj (2) 



From Eqnlll p = t^zjt . Substituting in Eqnpl 



The only unknown is n, and although this equation is hard to solve in closed form, the right hand side behaves 
linearly in the asymptotic sense. Using simple calculus techniques, we can estimate the solution to Eq|3]as: 

-2 
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The asymptotic approximation is so good that the rounding of n to an integer is the biggest source of error even 
for small values, say 10, of N]. Therefore, for practical purposes, Eq|4]gives the exact answer. To use this formula 
you have to know Nj, though. In spite of Ni being determined by N, it is not straightforward to find out an explicit 
formula. To obtain an approximation, we can simply generate a suitable number of networks using the BA algorithm 
and estimate Nr as the average of the number internal nodes over all the generated graphs. 

3 Simulations 

We programmed a simulator using about 700 lines of portable ANSI C. We ran a number of pseudo-random experi- 
ments. The main goal was to study the influence of the parameter a in the dynamics of the networks. More precisely, 
we were interested in answering, for a given network, the following questions as a function of a: 

1. Are there fixed points? 

2. If yes, 

(a) how many? 

(b) what fraction of trajectories converge to a fixed point? 

3. If we restrict the analysis to those trajectories that converged to a fixed point, 

(a) what is the average time (number of iterations) until a fixed point is reached? 

(b) what is the average work (total number of operations) [Juntil a fixed point is reached? 



Notice that, for a = 0, the amount of work per iteration is minimal (only 1 node is updated). On the other extreme, when a = 1, the amount 
of work per iteration is maximal (all the internal nodes are updated in parallel in each iteration). 



We setup the time horizon to 50000 iterations. This number seems to yield a robust determination of whether the 
network eventually reaches a fixed point or not. For the scale-free networks we ran the BA algorithm using parameters 
N = 100 and k = 2. We generated 50 networks. For each network we generated 1000 initial configurations. The states 
were generated randomly, but the values corresponding to source nodes were set to zero. More precisely, to generate 
the initial configuration s of a network TV, if K\ n = then the i-th component of s was set to zero. Otherwise, the i-th 
component of s was generated randomly by tossing a fair coin. 

Each one of the 50 x 1000 pairs network/initial configuration was used as a starting point for 1 1 simulations, each 
one using a different value of a. The values of a were 0, 0.1, ..., 0.9, 1, We ran each one of the 50 x 1000 x 11 
simulations until the trajectory converged to fixed point or the time horizon was reached. 

For the networks using the Erdos-Renyi topology we had to compute the input parameters for the ER algorithm 
(n and p). To estimate Nj, we generated 10000 graphs using the BA algorithm, with k = 2 and N = 100. Then 
we divided the total number of internal nodes by 10000, which yielded Nj = 62.745. By using EqHJ we determined 
n = 72. From Eq[T]we obtained p = 0.028196. 

With the two parameters, we repeated the process we used for the scale-free networks: We created 50 random 
networks and, for each one, we generated 1000 initial conditions. We also used the same values for a and the same 
time horizon for the simulations. We verified that the BA algorithm indeed produces graphs with vertex degrees that 
follow a power law probability distribution. Figurefllshows the relative frequency of each degree, for the 50 scale-free 
networks. The diagram, in log-log coordinates, shows a high similarity with a straight line. This is what we expect 
from a power-law distribution. The noise for the higher vertex degrees is not surprising, as highly connected vertices 
are rare and hence the population size is small. 
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Figure 1: Relative frequency of vertex degrees in the 50 scale-free networks 



4 Results and analysis 



4.1 Number of fixed points 

For both the scale-free and the Erdos-Renyi families we computed the number of fixed points found per each network 
and we summarize the results in Figures|2]and[3] The vertical axis of the histograms represent relative frequency over 
the 50 networks tested. Note that in 15 out of the 50 generated scale-free networks, no fixed point was found. More 
precisely, for all a, none of the 1000 trajectories was absorbed into a fixed point within the 50000 iterations used as 
time horizon. The number of Erdos-Renyi networks for which no fixed-point was found is 23. 

We proved analytically in Section|2]that the expected number of fixed points for any network is 1, Unfortunately, 
the actual distribution is hard to compute and therefore we use histograms as a base for a numerical approximation. 
We can see in Figure [2] that for scale-free network the most likely number of fixed points is 1. If we use the relative 





Figure 2: Frequency of number of fixed points for scale- Figure 3: Frequency of number of fixed points for Erdos- 
free networks. Renyi networks. 



frequencies to approximate probabilities, then the estimated expected number of fixed points is 1.08, which is close to 
the theoretical prediction. 

By contrast, for the Erdos-Renyi model, we can see in Figure [3] that the distribution gets more skewed and 
becomes the most likely value. The estimation of expected number of fixed points in this case yielded 0.88. This 
estimation is close to 1 if we consider we are using only 50 networks and the number of potential fixed points is about 
2 Nl , with N] sa 63, as we described in Subsection 2.2 



4.2 Fraction of converging trajectories 

We begin by focusing our analysis on those 35 scale-free networks for which we were able to prove the existence 
of fixed points (by finding them). Recall that we used 1000 initial configuration per network. If a < 1 (i.e, a — 
0, 0.1, . . . , 0.9) out of the 350000 stochastic trajectories we computed, 100% of them reached a fixed point. Despite 
the fact that for all these 35 networks we were also able to find at least one converging trajectory when a = 1, the 
situation in this case changed dramatically. In fact, out of the 35000 deterministic trajectories we simulated for those 
35 networks, only ~28.9% of them reached a fixed point. Notice also that the percentage of deterministic trajectories 
that reached a fixed point varied greatly form one network to another. We present the percentages for each network 
in Figure B] Every single trajectory that did not reach a fixed point ended in a cycle within the time horizon. For 
presentation purposes, we included the analogous diagram for a < 1 in Figure HI 

We immediately notice that a higher proportion of simulations reached a fixed point when a < 1 than when a = 1. 
This is not unexpected as the non-determinism in the stochastic trajectories is a way to avoid the cyclic attractors the 
deterministic trajectories tend to fall into. The remarkable property is that, if for a given network a single trajectory 
with a < 1 converged to a fixed point, then all the trajectories with a < 1 also converged (although not necessarily to 
the same point). Also note that all the fixed point discovered through fully parallel updates were also discovered with 
the a < 1 simulations. This indicates that stochastic updates can be used as a robust method to detect the existence of 
fixed points. 

For the Erdos-Renyi networks we computed the same statistics as we did for the scale-free networks, were ap- 
plicable. There were 27 out of the 50 networks were we found at least a fixed point. Interestingly, only 21 of those 
fixed points were detected by deterministic trajectories. The breakdown of percentages of converging trajectories by 
network for the case a = 1 is presented in Figure|7] Every single trajectory that did not reach a fixed point ended in a 
cycle within the time horizon. The breakdown of percentages by network for the case a < 1 is presented in Figure [6] 
We notice immediately the difference between Figures [6] and HI For some Erdos-Renyi networks with fixed point not 
all stochastic trajectories converged to a fixed point while all trajectories converged in the scale-free networks case, as 
we already described before. 

To emphasize the differences in the qualitative behavior we look more closely at two particular networks. We 
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Figure 4: Percentage of simulations that converged when 
a < 1 for each of the 35 scale-free networks with at least 
one fixed point. 



Figure 5: Percentage of simulations that converged when 
a — 1 for each of the 35 scale-free networks with at least 
one fixed point. 
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Figure 6: Percentage of simulations that converged for each 
of the 27 ER networks where a fixed point was found. 



Figure 7: Percentage of simulations that converged for each 
of the 2 1 ER networks where a fixed point was found. 



selected one Erdos-Renyi network for which the percentage of trajectories that converged (considering all a < 1) was 
about 50%. Recall that 1000 initial configurations were used, 10 values of a were tried and one trajectory per a was 
simulated. That means that about 5000 out of 10000 trajectories converged to a fixed point for that particular network. 
The histogram in Figure|9]shows that the choice of the initial configuration changes the probability of reaching a fixed 
point. The heights of the bars represent numbers of initial configurations. The horizontal axis is the percentage of 
trajectories (starting from one particular initial configuration) that converged to a fixed point. For instance, for the 
Erdos-Renyi network, we can deduce the following: 160 initial configurations converged in 40% of the simulations. 
Since for each initial configuration we ran exactly 10 simulations there is no need to associate intervals with bins. For 
contrast, we show in FigureJHlthe analogous histogram for an arbitrary scale-free network with a fixed point. 

Another difference between the influence of parallelism when updating scale-free or Erdos-Renyi networks is seen 
comparing Figures|2J[3]and[T0] Figure[l0]shows the difference in number of fixed points found depending on a being 1 
or less than 1. The columns that describe the a < 1 case represent the same values found in the histogram in Figure[3] 
The columns that describe the a = 1 case in Figure 10 show that fewer fixed points were discovered. This, again, 
strongly suggests that using partially parallel updates is a sensible strategy to find fixed points. It also indicates that 
the set of discovered fixed points somehow depends on the choice of a. For scale-free networks, we also discriminated 
between a = 1 and a < 1 cases. For these networks though, in both cases the histograms were identical to the one 
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Figure 8: Histogram of number of initial configurations Figure 9: Histogram of number of initial configurations 
that lead to a converging trajectory for a single scale-free that lead to a converging trajectory for a single Erdos-Renyi 
network which had at least one fixed point (a < 1). network which had at least one fixed point (a < 1). 

in Figure 13] This indicates that the convergence to fixed points is fairly insensitive to the choice of a. This property 
shows a form of robustness of scale-free networks. 
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Figure 10: Frequency of number of fixed points for ER networks for the a — 1 and a < 1 cases. 



4.3 Time and work until convergence 



For scale-free networks, Figures 11 and 12 represent the mean and maximum time/work to reach the fixed point for 
only the simulations that reached it. This means, for a < 1, each column gives the average and maximum over 35000 
simulations, while for a — 1 the values were computed over a smaller sample because only ^28.9% of the runs 
ended before the time horizon was reached. This makes the values of the rightmost column of each figure somewhat 
incomparable with the others. 



In Figure 1 1 we note how the average time changes with a. The qualitative behavior of the results are in part 



intuitive. If a = we expect the time to be high, since not much work is done per iteration. We also expected the time 
to increase when a approaches 1, because the system becomes more deterministic and tends to imitate the behavior 
of the a = 1 case. Therefore, the cycles of the deterministic trajectories become metastable regions of the stochastic 
trajectories. 

Note that the average number of iterations for a — 1 is smaller than for all the values for a < 1. The comparison 
is not fair, though, because the parallel update rule is not a reliable way to find fixed points. With respect to the amount 
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Figure 1 1 : Convergence time for scale-free networks 
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Figure 12: Convergence work for scale-free networks 
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Figure 13: Filtered convergence time for scale-free net- Figure 14: Filtered convergence work for scale-free net- 
works works 



of work, we can notice, as expected, that small values of a become very competitive (see Figure [12|). One possible 
explanation for the parallel simulations outperforming the a < 1 cases (neglecting the lack of robustness, of course) is 
as follows: The a = 1 updates are prone to solve only the "simple instances" of the problem. That is, only very stable 
fixed points located close to the initial configurations are likely to be found. 

To make a fair comparison between the a < 1 and a = 1 situations, we recompute the statistics but only for the 
cases where the deterministic model found a fixed point. This means, we considered only the pairs (network G, initial 



configuration s) such that the deterministic trajectory of G starting from s reached a fixed point. Figures 13 and 14 



show the results. From comparison against Figures 1 1 and 12 we notice that the average times of the simulations for 
a < 1 actually decreased. This supports the hypothesis of a = 1 working mostly for simple instances of the problem. 
The differences in running time may not be considered substantial though and the problem requires further study. 



For Erdos-Renyi networks, Figures 15 and 16 are the mean and maximum times/work to reach the fixed point for 
only the simulations that reached it. This means, for a < 1, each column gives the average and maximum over 27000 
simulations, while for a = 1 the values were computed over a smaller sample because only ^24% of the runs ended 
before the time horizon was reached. In this case we did not compute the averages/maximums for the stochastic model 
restricting the networks and initial conditions to those that reached a fixed point with deterministic updates. There was 



no qualitative difference between the results for scale-free networks (Figures 1 1 and 12 1 and the times for Erdos-Renyi 
networks. Data suggests that Erdos-Renyi networks are slower to converge than scale-free systems. 
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Figure 15: Convergence time for Erdos-Renyi networks Figure 16: Convergence work for Erdos-Renyi networks 

4.4 Discussion 

We will first make the case that the RBNs, as defined here, have a behavior that resembles that of natural networks. 
It is well known that gene expression/repression changes according to external cell conditions, such as temperature 
or availability of nutrients and/or oxygen. For brevity, we will refer to such conditions as EC, This can be modeled 
easily by setting the values of input nodes depending on the EC. It is also widely accepted that gene expression state 
is (after the transient) mostly function of the EC. In terms of the model this means the fixed point should be function 
of the values of the input nodes. Although stochasticity sometimes plays an important role in cell dynamics [25 1, for 
purposes of discussion we will assume that for each set of EC there should be one or a small number of fixed points. 
Another property is stability, meaning the fixed point should have a large basin of attraction. Additionally, the fixed 
point should be reached "quickly." By quickly, we mean that the number of transitions should be small compared to 
the total number of possible configurations. Finally, the behavior of the model should not depend dramatically on the 
choice of a. 

Our simulations with asynchronous updates on scale-free networks showed all properties mentioned above. The 
theorem in Subsection |2.1| together with the histogram in Figure [2] address the number of fixed points is, with high 
probability, close to one. The fact that all trajectories converged, as shown in Figure HI shows the model is (ro- 
bust/stable?) regardless of the choice of a, as long as it is less than 1. With respect to the speed of convergence, 
the average number of iterations to reach a fixed-point (see Figure) is small compared to the number of possible 
configurations, which was about 2 63 in all the numerical experiments. 

Of course, the goal of the model is to provide insight on biological process. The most remarkable observation is 
that the three aforementioned properties appeared in the simulations without requiring any kind of deliberate design. 
The network topologies, the dynamics and vertices to be updates were selected at random. Yet, the behavior was as 
desired, as long as the network was scale-free and a < 1. This suggests that the prevalent scale-free topologies in 
natural GRNs 1 3 1 could be a major factor in robustness of living cells. 

As an aside, finding fixed points of a Boolean network is equivalent to finding satisfying assignments to Boolean 
formulas. This is a well known problem in Computer Science, and is believed to be computationally hard if we consider 



the worst case running time. However, the results shown in Figures 12 14 and 16 which implies fast convergence 



suggest that the asynchronous update rule is an efficient heuristic to find fixed points or satisfying assignments. 

5 Conclusions 

We analyzed the tendency to reach a fixed point, the number of iterations and the work needed to converge to it for 
different families of Boolean networks and update policies. We summarize our results as follows: 

1. Using partial parallelism (i.e, a < 1), the likehood of the trajectory to reach a fixed point increases. This 
happens regardless of topology. 
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2. For scale-free networks the choice of a is not very important from the convergence point of view, as long as it 
is less than one. We base this conclusion on the fact that we could not find a simple example of a network not 
converging to a fixed point provided that it has 1 . 

3. For scale-free networks the choice of a is not very important from the work point of view, as long as it is less 
than one. This follows from the analysis of the filtered case. The work is about the same as in the a = 1 case, 
while full parallelism increases the likehood of being trapped in a cycle. 

4. For Erdos-Renyi networks the choice of a is relevant from the work point of view, when it is less than 1. The 
set of points discovered seems to be dependent somehow on the choice of a. 

5. The combination of scale-free networks and using a < 1 for the updates was enough to impose dynamical 
properties which are similar to those observed in nature. 
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